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The inclusive production rate of neutral pions in the rapidity range greater than y = 8.9 has been 
measured by the Large Hadron Collider forward (LHCf) experiment during y^s = 7TeV proton- 
proton collision operation in early 2010. This paper presents the transverse momentum spectra of 
the neutral pions. The spectra from two independent LHCf detectors are consistent with each other 
and serve as a cross-check of the data. The transverse momentum spectra are also compared with 
the predictions of several hadronic interaction models that are often used for high-energy particle 
physics and for modeling ultra-high-energy cosmic-ray showers. 

PACS numbers: 13.85.Tp, 13.85.-t 



I. INTRODUCTION 

One of the important tasks of strong-interaction 
physics described by Quantum Chromodynamics (QCD) 
is to provide a detailed understanding of forward particle 
production in hadronic interactions. QCD involves two 
types of limiting processes: "hard" and "soft" . 

Hard processes occur in the range characterized by 
a large four-momentum transfer i, where \t\ should be 
larger than 1 GcV 2 . Note that units used in this report 
are c = k (Boltzmann constant) = 1. Deep inelastic 
scattering that is accompanied by the exchange of vir- 
tual photons or vector bosons, or jets produced by large 
transverse momentum (j>t) partons are typical phenom- 
ena that are categorized as hard processes. The hard 
processes have been successfully described by perturba- 
tion theory, owing to the asymptotic freedom of QCD at 
high energy. 

On the other hand, soft processes occur when the four- 
momentum transfer \t\ is smaller than lGcV 2 . These 
processes, which correspond to a large impact parame- 
ter, have a large QCD coupling constant and cannot be 
calculated by perturbative QCD. Gribov-Regge theory 
is applicable for describing soft processes [l|, [2| and the 
Pomeron contribution, as a component of the Gribov- 
Regge approach to high-energy hadronic interactions, in- 
creases with increasing energy Q and should dominate 
at the TcV energy scale. However there still exists a 



problem for the theories that involve these virtual quasi- 
particles. Since the treatment of the Pomeron differs 
amongst the model theories they predict different results 
for particle production. Thus a deeper understanding 
of soft processes is needed and soft processes are mostly 
equivalent to forward or large rapidity particle produc- 
tion in hadronic interactions. However experimental data 
for large rapidity are meager. Moreover the experimental 
data that do exist have so far been carried out at rela- 
tively low energy, for example ISR [|| at ^/s = 53GeV 
and UA7 @ at ^fs = 630 GeV. 

The Large Hadron Collider forward (LHCf) experi- 
ment [(| has been designed to measure the hadronic pro- 
duction cross sections of neutral particles emitted in very 
forward angles in proton-proton collisions at the LHC, in- 
cluding zero degrees. The LHCf detectors have the capa- 
bility for precise measurements of forward high-energy 
inclusive-particle-production cross sections of photons, 
neutrons, and possibly other neutral mesons and baryons. 
Among the many secondary neutral particles that LHCf 
can detect, the tt° mesons are the most sensitive to the 
details of the proton-proton interactions. Thus a high 
priority has been given to analyzing forward ir° produc- 
tion data in order to provide key information for an as 
yet un-established hadronic interaction theory at the TeV 
energy scale. The analysis in this paper concentrates on 
obtaining the inclusive production rate for ir s in the ra- 
pidity range larger than y — 8.9 as a function of the w° 
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transverse momentum. 

In addition to the aim described above, this work is 
also motivated by an application to the understanding 
of Ultra-High-Encrgy Cosmic-Ray (UHECR) phenom- 
ena, which are sensitive to the details of soft ir° pro- 
duction at extreme energy. It is known that the lack of 
knowledge about forward particle production in hadronic 
collisions hinders the interpretation of observations of 
UHECR @, @. Although UHECR observations have 
made notable advances in the last few years crit- 
ical parts of the analysis depend on Monte Carlo (MC) 
simulations of air shower development that are sensitive 
to the choice of the hadronic interaction model. It should 
also be remarked that the LHC has reached 7TeV col- 
lision energy, which in the laboratory frame of UHECR 
observations is equivalent to 2.6 x 10 16 eV, and this en- 
ergy is above the "knee" region of the primary cosmic ray 
energy spectrum (~ 4x 10 15 cV) The data provided 
by LHCf should then provide a useful bench mark for the 
MC codes that are used for the simulation of UHECR at- 
mospheric showers. 

This paper is organized as follows. In Sec.|H]the LHCf 
detectors are described. Sec. IIIII summarizes the condi- 
tions for taking data and the MC simulation methodol- 
ogy. In Sec. HVI the analysis framework is described. The 
factors that contribute to the systematic uncertainty of 
the results are explained in Sec. [V] and the analysis re- 
sults are then presented in Scc. lVIl Sec. I VIII discusses the 
results that have been obtained and compare these with 
the predictions of several hadronic interaction models. 
Finally, concluding remarks are found in Sec. IVIIII 



II. THE LHCF DETECTORS 

Two independent LHCf detectors, called Arml and 
Arm2, have been installed in the instrumentation slots of 
the target neutral absorbers (TANs) [TtJ located ±140 m 
from the ATLAS interaction point (IP1) and at zero de- 
gree collision angle. Fig. [1] shows schematic views of the 
Arml (left) and Arm2 (right) detectors. Inside a TAN 
the beam-vacuum-chamber makes a Y-shaped transition 
from a single common beam tube facing IP1 to two sepa- 
rate beam tubes joining to the arcs of the LHC. Charged 
particles produced at IP1 and directed towards the TAN 
are swept aside by the inner beam separation dipole mag- 
net Dl before reaching the TAN. Consequently only neu- 
tral particles produced at IP1 enter the LHCf detector. 
At this location the LHCf detectors cover the pseudo- 
rapidity range from 8.7 to infinity for zero degree beam 
crossing angle. With a maximum beam crossing angle of 
140 yurad, the pseudorapidity range can be extended to 
8.4 to infinity. 

Each LHCf detector has two sampling and imaging 
calorimeters composed of 44 radiation lengths (Xo) of 
tungsten and 16 sampling layers of 3 mm thick plastic 
scintillator. The transverse sizes of the calorimeters are 
20 x20mm 2 and 40 x40mm 2 in Arml, and 25 x25mm 2 



and 32 x 32 mm 2 in Arm2. The smaller calorimeters cover 
zero degree collision angle. Four X-Y layers of position 
sensitive detectors are interleaved with the layers of tung- 
sten and scintillator in order to provide the transverse 
positions of the showers. Scintillating fiber (SciFi) belts 
are used for the Arml position sensitive layers and silicon 
micro-strip sensors are used for Arm2. Readout pitches 
arc 1mm and 0.16 mm for Arml and Arm2, respectively. 

More detail on the scientific goals and the construc- 
tion and performance of the detectors can be found in 
previous reports [l8|-[22j. 




FIG. 1: (color online). Schematic views of the Arml (left) 
and Arm2 (right) detectors. The transverse sizes of the 
calorimeters are 20 x20mm 2 and 40 x40mm 2 in Arml, and 
25 x25mm 2 and 32 x32mm 2 in Arm2. 



III. SUMMARY OF THE CONDITIONS FOR 
TAKING DATA AND OF THE METHODOLOGY 
FOR PERFORMING MONTE CARLO 
SIMULATIONS 

A. Conditions for taking experimental data 

The experimental data used for the analysis of this 
paper were obtained on May 15th and 16th 2010 during 
proton- proton collisions at y / s=7TeV with zero degree 
beam crossing angle (LHC Fill 1104). Data taking was 
carried out in two different runs: the first run was on May 
15th from 17:45 to 21:23, and the second run was on May 
16th from 00:47 to 14:05. The events that were recorded 
during a luminosity optimization scan and a calibration 
run were removed from the data set for this analysis. 

The range of total luminosity of the three crossing 
bunch pairs was C = (6.3 — 6.5) x 10 28 cm~ 2 s _1 for the first 
run and C = (4.8 — 5.9) x 10 28 em _2 s _1 for the second run. 
These ranges of luminosity were ideal for the LHCf data 
acquisition system. The integrated luminosities for the 
data analysis reported in this paper were derived from the 
counting rate of the LHCf Front Counters [23| , and were 
2.53 nb" 1 (Arml) and 1.90 nb" 1 (Arm2) after taking the 
live time percentages into account. The average live time 
percentages for the first/second run were 85.7%/81.1% 
for Arml and 67.0%/59.7% for Arm2. The live time 
percentages for the second run were smaller than the first 
run owing to a difference in the trigger schemes. In both 
runs the trigger efficiency achieved was >99% for photons 
with energy E > 100 GeV [24|. 
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The events containing more than one collision in a 
single bunch crossing (pile-up events) could potentially 
cause a bias in the pt spectra. For example combinatorial 
single-hits from different collisions within a single bunch 
crossing might be identified as multi-hit events from a sin- 
gle collision and removed from the analysis. (Multi-hit 
events have two showers in a single calorimeter and are 
eliminated from the data analysis. The production rates 
are later corrected for this cut. See Fig. [2] and related 
discussion.) However it can be shown that pile- up events 
are negligible for the LHCf data taking conditions of this 
report. Given that a collision has occurred, the proba- 
bility of pile-up (Ppiieup) is calculated from the Poisson 
probability distribution for n collisions P po i(^) according 
to Fpiicup = P poi (n > 2)/P poi (n > 1). With the high- 
est bunch luminosity C = 2.3 x 10 28 cm~ 2 s _1 used in this 
analysis, an inelastic cross section cri no i = 73.6 mb and the 
revolution frequency of LHC / rev = 11.2 kHz, the pile-up 
probability is P p ii CU p ~ 0.07. However considering that 
the acceptance of the LHCf calorimeter for inelastic col- 
lisions is ~0.03, only 0.2% of events have more than one 
shower event in a single calorimeter due to pile-up and 
this is negligible. 

Detailed discussions of pile-up effects and background 
events from collisions between the beam and residual gas 
molecules in the beam tube can be found in previous 
reports [g|,[24[. 



B. Methodology for performing Monte Carlo 
simulations 

MC simulation consists of three steps: (1) proton- 
proton interaction event generation at IP1, (2) transport 
from IP1 to the LHCf detectors and (3) the response of 
the LHCf detectors. 

Proton-proton interaction events at ^fs = 7 TeV and 
the resulting flux of secondary particles and their kine- 
matics are simulated with COSMOS (version 8.81). COS- 
MOS acts as the front end for the external hadronic in- 
teraction models (qgsjet 11-03 [25Tl. dpmjet 3.04 [26j |. 
SIBYLL 2.1 [13 and EPOS 1.99 [28j]) that describe the 
proton-proton interactions. While pythia 8.145 [29l.l3"bT| 
serves as its own front end for the generation of proton- 
proton interaction events. 

Next, the generated secondary particles arc trans- 
ported in the beam pipe from IP1 to the TAN, taking 
account of the deflection of charged particles by the Ql 
quadrupole and Dl beam separation dipolc, particle de- 
cay, and particle interaction with the beam pipe and the 
Y-shaped beam- vacuum-chamber transition made of cop- 
per (IXq projected thickness in front of the LHCf de- 
tectors). Charged particles are swept away by the Dl 
magnet before reaching the LHCf detectors. This sim- 
ulation uses the epics library [3l[ (version 7.49) and a 
part of COSMOS, epics deals with the transport of sec- 
ondary particles. Particle interactions with the residual 
gas molecules inside the beam pipe are not simulated. 



Contamination from beam-gas background events in the 
data set used for analysis is estimated to be only ~ 0.1 % 
and has no significant impact on the pt spectra reported. 

Finally the simulations of the showers produced in 
the LHCf detectors and their response arc carried out 
for the particles arriving at the TAN using the COSMOS 
and EPICS libraries. The survey data for detector po- 
sition and random fluctuations equivalent to electrical 
noise are taken into account in this step. The Landau- 
Pomeranchuk-Migdal effect [HI, [33| that longitudinally 
lengthens an electromagnetic shower at high energy is 
also considered. A change of the pt spectra caused by 
LPM effects is only at the 1% level since the reconstruc- 
tion of energy deposited in the calorimeters is carried out 
to a sufficiently deep layer whereby the energy of electro- 
magnetic showers is almost perfectly deposited within the 
calorimeter. 

The simulations of the LHCf detectors are tuned to 
test beam data taken at the CERN SPS in 2007 [1|. 
The validity of the detector simulation was checked by 
comparing the shower development and deposited energy 
for each calorimeter layer to the results obtained by the 
fluka library [34] ]. 

In order to validate the reconstruction algorithms and 
to estimate a possible reconstruction bias beyond the en- 
ergy range of the SPS test beam results, the MC sim- 
ulations are generated for 1.0 x 10 8 inelastic collisions, 
where the secondary particles are generated by the EPOS 
1.99 [28| hadronic interaction model. This MC simula- 
tion is referred to as the "reference MC simulation" in 
the following text. 

Similarly the "toy MC simulations" discussed below 
are performed in order to determine various correction 
factors to use in the event reconstruction processes. In 
the toy MC simulations, a single-photon with a given 
fixed energy is directed at the LHCf detectors. 



IV. ANALYSIS FRAMEWORK 

A. Event reconstruction and selection 

Observation of tt° mesons by a LHCf detector is illus- 
trated in Fig. [21 The ir°s are identified by their decay 
into two photons. Since the 7r°s decay very close to their 
point of creation at IP1, the opening angle (8) between 
the two photons is the transverse distance between pho- 
ton impact points at the LHCf detectors divided by the 
distance from IP1 (z ~ ±141.05 m). Consequently the 
opening angle for the photons from 7r° decay that are de- 
tected by a LHCf detector is constrained by 9 < 0.4 mrad 
for Arml and 9 < 0.6 mrad for Arm2. Other kinematic 
variables of the ir°s (energy, pp, and rapidity) are also 
reconstructed by using the photon energy and incident 
position measured by each calorimeter. Note that for the 
analysis of this paper events having two photons enter- 
ing the same calorimeter (multi-hit events) are rejected 
(right panel of Fig. [5]). The accuracy of energy recon- 



4 



struction for such events is still under investigation. The 
final inclusive production rates reported in this paper are 
corrected for this cut. In order to ensure good event re- 
construction efficiency, the range of the 7r° rapidity and 
Pt are limited to 8.9 < y < 11.0 and < 0.6 GeV, 
respectively. All particles other than photons from 7r° 
decay arc ignored in this analysis. Thus, also accord- 
ing to the multi-hit 7r° correction described in detail in 
Sec. IIVF1 the reported production rates are inclusive. 
The standard reconstruction algorithms are described in 
this section and systematic uncertainties will be discussed 
in Sec. El 




FIG. 2: (color online). Observation of ty° decay by a LHCf 
detector. (Left) Two photons enter different calorimeters. 
(Right) Two photons enter one calorimeter. 



1. Hit position reconstruction 

The transverse impact positions of particles entering 
the calorimeters are determined using the information 
provided by the position sensitive layers. In this analysis, 
the transverse impact position of the core of an electro- 
magnetic shower is taken from the position of the peak 
signal on the position sensitive layer that has the largest 
energy deposited amongst all the position sensitive lay- 
ers. 

Hit positions that fall within 2 mm of the edges of the 
calorimeters are removed from analysis due to the large 
uncertainty in the energy determination of such events 
owing to shower leakage. For the toy MC simulations, the 
position reconstruction resolution is defined as the one 
standard deviation difference between the true primary 
photon position and the reconstructed position of the 
shower axis. The estimated resolution using the toy MC 
simulations and test beam data for a single photon with 
energy E > 100 GcV is better than 200 am and 100 /im 
for Arml and Arm2, respectively [2ll.[35j. 

Multi-hit events defined to have more than one photon 
registered in a single calorimeter are eliminated from the 
analysis in this paper. Multi-hit candidates that have 
two distinct peaks in the lateral shower impact distribu- 
tion are searched for using the algorithm that has been 
implemented in the TSpectrum [3a] class in ROOT [37j . 
When the separation between peaks is greater than 1 mm 
and the lower energy photon has more than 5 % of the en- 



ergy of the nearby photon, the MC simulation estimated 
efficiencies for identifying multi-hit events are larger than 
70 % and 90 % for Arml and Arm2, respectively Q- The 
efficiency for Arm2 is better than that for Arml owing 
to the finer readout pitches of the silicon micro-strip sen- 
sors. The subtraction of the remaining contamination by 
multi-hit events is discussed in Sec. IIV CI 

On the other hand for single-hit events not having 
two identifiable peaks, the MC simulation estimated effi- 
ciency for correctly identifying true single photon events 
with energy E > 100 GeV is better than 98 % both for 
Arml and Arm2, although the precise percentage de- 
pends slightly on the photon energy. 



2. Energy reconstruction 

The charge information in each scintillation layer is 
converted to a deposited energy by using calibration fac- 
tors obtained from the SPS electron test beam data taken 
below 200 GeV pj| . In this analysis the deposited en- 
ergy is scaled to the number of minimum ionizing shower 
particles with a coefficient 1MIP = 0.453 MeV that cor- 
responds to the most probable deposited energy by a 
150 GeV muon passing through a 3 mm thick plastic scin- 
tillator. The sum of the energy deposited in the 2 nd to 
13 th scintillation layers (dE [MIP]) is then converted to 
the primary photon energy i?[GeV] using a polynomial 
function 

E = A E dE 2 + B E dE + C E - (1) 

The coefficients Ae [GeV/MIP 2 ], B E [GeV/MIP] and 
C E [GeV] are determined from the response of the 
calorimeters to single photons by the toy MC simula- 
tions. The validity of this method has been confirmed 
with the SPS beam tests. The MC estimated energy res- 
olution for single photons above 100 GeV considering the 
LHC data taking situation is given by the expression 

a(E)/E ~ 8 %/y/E/100 GeV © 1 %. (2) 

Corrections for shower leakage effects [TH, [l|| are car- 
ried out during the energy reconstruction process. Cor- 
rections are applied for leakage of particles out of the 
calorimeters and for leakage of particles in that have es- 
caped from the adjacent calorimeter. Both of the leakage 
effects depend on the transverse location of the shower 
axis in the calorimeters. The correction factors have been 
estimated from the toy MC simulations. The light-yield 
collection efficiency of the plastic scintillation layers [l9| 
is also a function of the transverse location of the shower 
axis and corrected for in this step. 

Events having a reconstructed energy below 100 GeV 
are eliminated from the analysis firstly to reject particles 
produced by interaction of collision products with the 
beam pipe, and secondly to avoid errors due to trigger 
inefficiency (see Sec. IIII A"|) . 
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3. Particle identification 

The particle identification (PID) process is applied in 
order to select pure electromagnetic showers, specifically 
photons from 7r° decay, and to reduce hadron contami- 
nation, specifically from neutrons. A parameter Lg % is 
defined for this purpose. L go % is the longitudinal dis- 
tance, in units of radiation length, measured from the 
1st tungsten layer of a calorimeter to the position where 
the energy deposition integral reaches 90 % of the total 
shower energy deposition. Fig. [3] shows the distribution 
of £90% for the 20 mm calorimeter of the Arml detec- 
tor for events having a reconstructed energy in the range 
500GcV< E <1TcV. Experimental data (black dots) 
and the MC simulations based on QGSJET 11-03 (shaded 
areas) are shown. The normalization factors of pure pho- 
ton and pure hadron incident events are modified to get 
the best agreement between the £90% distributions of the 
experimental data and the MC simulations. The best 
agreement is obtained by a chi-squarc test of the £90% 
distribution of the experimental data relative to the MC 
simulation. The two distinct peaks correspond to photon 
(£90% iS 20 Xq) and hadron (L 90 % > 20 X ) events. 

PID criteria that depend on the energy of the individ- 
ual photons are defined in terms of the £90% distribution 
in order to keep the ir° selection efficiency at approx- 
imately 90% over the entire px range. These criteria 
fi,90%(Ei, E2) are expressed as a function of the pho- 
ton energies measured by the small (Ei) and large (E2) 
calorimeters and have been determined by the toy MC 
simulations for each Arm. The remaining hadron con- 
tamination is removed by background subtraction intro- 
duced in Sec. IIVC] The unavoidable selection inefficiency 
of 10 % is corrected for in the unfolding process to be dis- 
cussed later (Sec. HVDj) . 

Table U summarizes the ir° event selection criteria that 
are applied prior to reconstruction of the 7r° kinematics. 




5 10 15 20 25 30 35 40 45 
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FIG. 3: (color online). L 90 % distribution measured by 
the Arml-20mm calorimeter for the reconstructed energy of 
500 GeV-1 TeV. 



Incident position 


within 2 mm from the edge of calorimeter 


Energy threshold 


^photon > 100 GeV 


Number of hits 


Single-hit in each calorimeter 


PID 


Photon like {L 90% < f L90 %{Ei, E 2 )) 



TABLE I: Summary of criteria for event selections of the 7r 
sample. 



B. ty° reconstruction 

Candidates for ir° events are selected using the charac- 
teristic peak in the two-photon invariant mass distribu- 
tion corresponding to the ir° rest mass. Reconstruction 
of the invariant mass m 77 is done using the incident po- 
sitions and energies information of the photon pair, 

™ 77 = (<7i+<7 2 ) 2 -EM 2 , (3) 

where qi and Ei are the energy-momentum 4- vectors and 
energies of the decay photons in the laboratory frame, 
respectively. 9 is the opening angle between the two 
photons in the laboratory frame. The last approxima- 
tion in Eq. (|3|) is valid since the 7T°s decay very close 
to IP1 (mean ir° flight path < 1mm). This approxi- 
mation and the reconstruction algorithm for 7r° events 
have been verified by analysis of the reference MC sim- 
ulations of the energy, rapidity and pt of the 7r°s. The 
reconstructed invariant mass is concentrated near peaks 
at 135.2±0.2MeV in Arml and 134.8±0.2MeV in Arm2, 
thus reproducing the ir° mass. The uncertainties given 
for the mass peaks are statistical only. 

It should be noted however that in the ir° analysis of 
the experimental LHCf data energy scale corrections are 
needed so the ir° mass peaks for Arml and Arm2 occur 
at the proper value. With no energy scale corrections 
applied to the LHCf data, the reconstructed invariant 
mass peaks using gain calibration constants determined 
by test beam data occur at 145.8±0.1 MeV (Arml) and 
139.9±0.1 MeV (Arm2). Therefore energy scale correc- 
tions of -8.1% (Arml) and -3.8% (Arm2) applied to 
the raw measured photon energies are needed to bring 
the reconstructed tt° rest mass into agreement with the 
world averaged it rest mass [38| . The cause of these en- 
ergy scale corrections is probably due to a temperature 
dependent shift of PMT gain. However at this point the 
temperature dependent shift of PMT gain is only qual- 
itatively understood. Note that the typical uncertainty 
in opening angle is estimated to be less than 1 % relative 
to the reconstructed invariant mass by the position de- 
termination resolution and the alignment of the position 
sensitive detectors. 



C. Background subtraction 

Background contamination of two-photon ir° events by 
hadron events and the accidental coincidence of two pho- 
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tons not coming from the decay of a single ir° are sub- 
tracted using the so called "sideband" method. 

Fig. H] shows an example of the reconstructed two- 
photon invariant mass distribution of the experimental 
data of Arml in the rapidity range from 9.0 to 9.2. The 
energy scale correction discussed in the previous section 
has been applied. The sharp peak around 135 McV is due 
to 7T° events. The solid curve represents the best-fit of a 
composite physics model to the invariant mass distribu- 
tion of the data. The model consists of an asymmetric 
Gaussian distribution (also known as a bifurcated Gaus- 
sian distribution) for the signal component and a 3rd or- 
der Chebyshev polynomial function for the background 
component. The dashed curve indicates the background 
component. 

Using the expected mean {rh) and 1 a deviations (07 for 
lower side and a u for upper side) of the signal component, 
the signal window is defined as the invariant mass region 
within the two solid arrows shown in Fig. 21 where the 
lower and upper limits are given by m — 3cr; and m + 
3<t u , respectively. The background window is constructed 
from the two sideband regions, [rh — 6(7;, rh — 3a{\ and 
[m + 3cr u , TO + 6er u ], that arc defined as the invariant mass 
regions within the dashed arrows in Fig. [4] 

The rapidity and pt distributions of the signal 
(/(y,Pr) ) are then obtained by subtracting the back- 
ground distribution (/(y,pr) BG ), estimated by the 
background window, from the signal-rich distribution 
(/(ZA Pt) Sis+BG ) selected from the signal window. The 
fraction of the background component included in the 
signal window can be estimated using the likelihood func- 
tion (Lbg(j/ 5 Pt, m 'yy)) characterized by the best-fit 3rd 
order Chebyshev polynomial function. For simplicity, 
LBG(y,PTi m ~fy) is shortened as Lbg in the following 
text. Thus the signal distribution with background sub- 
tracted is given by 

/(y,p T ) Sig = f(y, P ?) S[g+BG - (4) 

R(y,PT, rh, ai,(T u )f(y,p T ) BG , 

where R(y,pT,rh,ai,a u ) is the normalization for the 
background distribution and written as 

R(y,PT,m,ai,a u ) = (5) 

Jm-3 CTi L BG^77 

f m-3(T, T , r m+6a u T , 

Jfh-Gv, L BGdWi 7 7 + Jrh+3a u L BGdm 71 



D. Unfolding of spectra 

The raw rapidity - px distributions must be corrected 
for unavoidable reconstruction inefficiency and for the 
smearing caused by finite position and energy resolutions. 
An iterative Bayesian method [H, H(| is used to simul- 
taneously correct for both effects. The advantages of an 
iterative Bayesian method with respect to other unfold- 
ing algorithms are discussed in another report [39| . The 
unfolding procedure for the data is organized as follows. 




Reconstructed [MeV] 

FIG. 4: (color online). Reconstructed invariant mass distri- 
bution within the rapidity range from 9.0 to 9.2. Solid curve 
shows the best-fit composite physics model to the invariant 
mass distribution. Dashed curve indicates the background 
component. Solid and dashed curves indicate the signal and 
background windows, respectively. 



First, the response of the LHCf detectors to single 7r° 
events is simulated by toy MC calculations. In the toy 
MC simulations, two photons from the decay of 7r°s and 
low energy background particles such as those originat- 
ing in a prompt photon event or a beam-pipe interac- 
tion are traced through the detector and then recon- 
structed with the event reconstruction algorithm intro- 
duced above. Note that the single ir° kinematics that 
are simulated within the allowed phase space are inde- 
pendent of the particular interaction model that is be- 
ing used. The background particles are simulated by a 
hadronic interaction model which is discussed later, since 
the amount of background particles is not directly mea- 
sured by the LHCf detector. 

The detector response to 7T° events depends on rapidity 
and pt, since the performance of the particle identifica- 
tion algorithm and the selection efficiency of events with 
a single photon hit in both calorimeters depend upon the 
energy and the incident position of a particle. The recon- 
structed rapidity - pt distributions for given true rapid- 
ity - pt distributions then lead to the calculation of the 
response function. Then the reconstructed rapidity and 
Pt spectra are corrected with the response function which 
is equivalent to the likelihood function in Bayes' theo- 
rem. The corrections are carried out iterativcly whereby 
the starting point of the current iteration is the ending 
point of the previous iteration. Statistical uncertainty is 
also propagated from the first iteration to the last. Iter- 
ation is stopped at or before the 4th iteration to obtain 
a regularization of the unfolded events. 

Validation of the unfolding procedure is checked by 
applying the response function to the reference MC sim- 
ulation samples. The default response function is deter- 
mined with two photons from ir° decay and the low en- 
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ergy (E < 100 GcV) background particles generated by 
epos 1.99. Validity of the choice of epos 1.99 is tested 
by comparing two corrected spectra, one generated by 
EPOS 1.99 and another by pythia 8.145. No statisti- 
cally significant difference between the corrected spectra 
is found. A chi-square test of the corrected spectra based 
on the default response function against the true spectra 
ensures the chi-square probability is greater than 60%. 
Thus it is concluded that with the background subtrac- 
tion and unfolding methods used in this analysis there 
is no significant bias and the statistical uncertainty is 
correctly quoted. Accordingly no systematic uncertainty 
related to the choice of the hadronic interaction mod- 
els for the reference MC simulations is considered in the 
analysis that follows. 



E. Acceptance and branching ratio correction 

The apertures of the LHCf calorimeters do not cover 
the full 2tt azimuthal angle over the entire rapidity range 
that is sampled. A correction for this is applied to the 
data before it is compared with theoretical expectations. 

The correction is done using the rapidity - phase 
space. Correction coefficients are determined as follows. 
First, using a toy MC simulation, a single 7r° is generated 
at IP1 and the decay photons are propagated to the LHCf 
detectors. The energy-momentum 4-vectors of the 7r°s 
are randomly chosen so that they cover the rapidity range 
that the LHCf detectors are able to measure. The beam 
pipe shadow on the calorimeter and the actual detector 
positions are taken into account using survey data. 

Next fiducial area cuts in the transverse X-Y plane 
are applied to eliminate particles that do not fall within 
the acceptance of the calorimeters. In the fiducial area 
cuts, a systematic shift of the proton beam axis is applied 
according to the reconstruction of the beam-axis during 
LHC operation. In addition a cut is applied to eliminate 
photons with energy less than 100 GeV. This corresponds 
to the treatment of the actual data for reducing the back- 
ground contamination by particle interactions with the 
beam pipe. 

Finally two phase space distributions of 7r°s are pro- 
duced; one is for all 7r°s generated at IP1 and the other 
is for 7r°s accepted by the calorimeters. The ratio of the 
distribution of accepted 7r°s divided by the distribution 
of all 7r°s is then the geometrical acceptance efficiency. 
Fig. [5] shows the acceptance efficiency as a function of 
the 7T° rapidity and pt and dashed curves indicate lines 
of constant tt° energy E = 1 TcV, 2 TeV and 3 TcV. The 
left and right panels indicate the acceptance efficiency 
for Arml and Arm2, respectively. The final rapidity and 
Pt spectra are obtained by applying the acceptance map 
shown in Fig.[S]to the acceptance uncorrected data. Note 
that the correction maps in Fig. [5] are purely kinematic 
and do not depend upon the particular hadronic interac- 
tion model that has been used. The uncertainty of the 
acceptance map caused by the finite statistics of the MC 



simulations is negligible. 

The branching ratio of ir° decay into two photons is 
98.8% and then inefficiency due to 7r° decay into chan- 
nels other than two photons (1.2 %) is taken into account 
by increasing the acceptance efficiency in rapidity - pt 
phase space by 1.2% everywhere and is independent of 
the particular hadronic interaction model. 




Rapidity 



Rapidity 



FIG. 5: (color online) . The acceptance map of tv° detection by 
the LHCf detectors in rapidity - px phase space: Arml (left) 
and Arm2 (right). Fiducial area cuts and energy threshold 
(-Ephoton > 100 GeV) are taken into account. Dashed curves 
indicate lines of constant energy 7r°s, E = 1 TeV, 2 TeV and 
3 TeV. 



F. Multi-hit 7T° correction 

The detected events have been classified into two types 
of events: single-hit ir° and multi-hit 7r° events. The for- 
mer class consists of two photons, one in each of the 
calorimeters of an Arml or Arm2 detector. A multi- 
hit 7T° event is defined as a single ir° accompanied with 
at least one additional background particle (photon or 
neutron) in one of the calorimeters. In this analysis, 
only single-hit ir° events are considered, and multi-hit 
7r events are rejected in the single-hit selection process 
(Sec. IIV A 1[) when the energy of the additional back- 
ground particle is beyond the energy threshold of the 
cut. 

The loss of multi-hit ir° events is corrected for with the 
help of event generators. A range of ratios of multi-hit 
plus single-hit to single-hit ir° events is estimated us- 
ing several hadronic interaction models in each rapidity 
range. The observed pt spectra are then multiplied by 
the mean of these ratios and also contribute a systematic 
uncertainty corresponding to the variation among the in- 
teraction models. In this way the single-hit tt° spectra 
are corrected so they represent inclusive ir° production 
spectra. The pt dependent range of the flux of multi-hit 
7T° events has been estimated using QGSJET 11-03, DPM- 
jet 3.04, SIBYLL 2.1, EPOS 1.99 and pythia 8.145, and 
resulted in a range of %-10 % of the flux of single- hit 
7T° events. 
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V. SYSTEMATIC UNCERTAINTIES 

A. Energy scale 

The known rest mass of the 7r°s is 134.9766 ± 
0.0006 MeV [38[ whereas the peak of the two-photon in- 
variant mass measured by the two LHCf detectors occurs 
at 145.8±0.1McV (Arml) and 139.9±0.1MeV (Arm2) 
where the ± 0.1 MeV uncertainties are statistical. The 
mass excess error is +8.1% for Arml and +3.8% for 
Arm2. According to Eq. [3] there are two possible sources 
for mass excess error; (1) systematic over estimates of 
the energies E\ and E2 of the two decay photons and (2) 
systematic over estimate of the opening angle between 
the two photons. As discussed in Sec. IV B the typical 
uncertainty in opening angle is less than 1 %, too small 
to explain the observed mass excesses. This leaves mea- 
surement of the photon energies as the source of mass 
excess error. 

The uncertainty in measurement of photon energy has 
also been investigated in a beam test at SPS and calibra- 
tion with a radiation source. The estimated uncertainty 
of photon energy from these tests is 3.5 %. The 3.5 % un- 
certainty is dominated by the uncertainties in factors con- 
verting measured charge to deposited energy [20j . Note 
that the linearity of each PMT was carefully tested be- 
fore detector assembly over a wide range of signal am- 
plitude by exciting the scintillator with a 337 nm UV 
laser pulse [f| [l!| . The difference of reconstructed energy 
between the reconstruction algorithm with and without 
non-linearity correction of PMTs for 3TeV photons is 
only 0.5% at maximum, nevertheless the measured non- 
linear response functions have been applied in the anal- 
ysis. 

The systematic uncertainties estimated by the beam 
test data at SPS (3.5 % for both Arms) are considered 
as uncorrelated among the pt bins, while the systematic 
uncertainties owing to the mass excess errors (8.1% for 
Arml and 3.8 % for Arm2) are considered as correlated 
between each pt bin. The systematic shift of bin contents 
due to the energy scale uncertainties is estimated using 
two energy spectra by artificially scaling the energy with 
the two extremes. The ratios of the two extreme spectra 
to the non-scaled spectrum are assigned as systematic 
shifts in each bin. 



B. Particle identification 

The L 90 % distribution described in Sec. IIV A~3l is used 
to select LHCf ir° events for the pt spectra presented in 
Sec. lVIl Some disagreements in the L 90 % distribution are 
found between the LHCf data and the MC simulations. 
This may be caused by residual errors of the channel-to- 
channcl calibrations of the LHCf detector relative to the 
LHCf detector simulation. 

The corresponding systematic uncertainty of the £90% 
distribution is evaluated by comparing the Lg % distri- 



bution of the LHCf 7T° candidate events of the measured 
data with the MC simulation. The L 90 o/ distribution for 
LHCf 7T° events is increased by at most one radiation 
length compared to the MC simulation. The systematic 
shifts of pt spectra bin contents are taken from the ratio 
of Pt spectra with artificial shifts of the Lg % distribution 
to the pt spectra without any L 90 % shift. This effect may 
distort the measured pt spectra by 0-20 % depending on 
Pt- 

C. Offset of beam axis 

In the geometrical analysis of the data, the projected 
position of the zero degree collision angle at the LHCf 
detectors (beam center) can vary from fill to fill owing to 
slightly different beam transverse position and crossing 
angles at IP1. The beam center at the LHCf detectors 
can be determined by two methods; first by using the 
distribution of particle incident positions measured by 
the LHCf detectors and second by using the information 
from the Beam Position Monitors (BPMSW) installed 
±21 m from IP1 [4l|. Consistent results for the beam 
center are obtained by the two methods applied to LHC 
fills 1089-1134 within 1mm accuracy. The systematic 
shifts to pt spectra bin contents are evaluated by taking 
the ratio of spectra with the beam-center displaced by 
1 mm to spectra with no displacement as determined by 
the distribution of particle incident positions measured 
by the LHCf detectors. Owing to the fluctuations of the 
beam-center position, the pt spectra are modified by 5— 
20 % depending on the rapidity range. 

D. Single- hit selection 

Since energy reconstruction is degraded when more 
than one photon hits a given calorimeter, only single-hit 
events are used in the analysis. Owing to selection effi- 
ciency greater than 98 % for single-hit events and rejec- 
tion of contamination by multi-hit events by the invariant 
mass cut, the systematic shift caused by the uncertainty 
in single-hit selection to bin contents is 3%. 

E. Position dependent correction 

As described in Sec. IIV A2[ energy reconstruction of 
the photons is sensitive to shower leakage effects which 
are a function of the photon incident position. System- 
atic uncertainties related to the leakage-out and leakagc- 
in effects arise from residual errors of calorimeter re- 
sponse when tuning of the LHCf detector simulation to 
the calibration data taken at SPS [20| that then lead to 
a mis-reconstruction of energy. Another source of uncer- 
tainties in energy reconstruction is an error in light-yield 
collection efficiency which is also dependent on the pho- 
ton incident position. 
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The systematic uncertainty due to position dependent 
effects is estimated by comparing two distributions of the 
energy deposited at each incident position bin. The first 
distribution is taken from the beam tests at SPS and the 
second distribution is generated by toy MC simulations 
that assume the upstream geometry of the test beam at 
SPS. Shifts of reconstructed px attributed to the residual 
errors in calorimeter response between these two energy 
distributions are assigned as the systematic uncertainties. 
The typical systematic shifts of Arml (Arm2) are 5 % 
(5 %) for low p T and 40 % (30 %) for large px- Owing to 
the light guide geometry, the systematic uncertainty of 
the Arml detector is larger than the Arm2 detector. 



F. Luminosity 

The instantaneous luminosity is derived from the 
counting rate of the Front Counters (FC). The calibration 
of the FC counting rates to the instantaneous luminosity 
was made during the Van der Meer scans on April 26th 
and May 9th 2010 [23[. The calibration factors obtained 
from two Van der Meer scans differ by 2.1 %. The esti- 
mated luminosities by the two FCs for the May 15th data 
differ by 2.7%. Considering the uncertainty of ±5.0% in 
the beam intensity measurement during the Van der Meer 
scans [lU, we estimate an uncertainty of ±6.1% in the 
luminosity determination. 



VI. RESULTS OF ANALYSIS 

The px spectra derived from the independent analyses 
of the Arml and Arm2 detectors are presented in Fig. m 
for six ranges of rapidity y: 8.9 to 9.0, 9.0 to 9.2, 9.2 to 
9.4, 9.4 to 9.6, 9.6 to 10.0 and 10.0 to 11.0. The spectra 
in Fig. [6] are after all corrections discussed in previous 
sections have been applied. The inclusive production rate 
of neutral pions is given by the expression 

J—E— = 1 d 2 N( PT ,y) 

dp 3 N inei 2ir ■ p T ■ dp T ■ dy ' 

<7i ne i is the inelastic cross section for proton-proton colli- 
sions at y/s — 7TeV. Ed 3 a /dp 3 is the inclusive cross sec- 
tion of 7T° production. The number of inelastic collisions, 
iVjnei, used for normalizing the production rates of Fig. [5] 
has been calculated from A^ ne i = <?mei / £dt, assuming 
the inelastic cross section a- lILe \ = 73.6 mb. This value for 
(Ti no i has been derived from the best COMPETE fits [H| 
and the TOTEM result for the elastic scattering cross 
section [43|. Using the integrated luminosities reported 
in Sec. |IIT3 N inei is 1.85xl0 8 for Arml and 1.40xl0 8 
for Arm2. d 2 N(pT,y) is the number of 7r°s detected in 
the transverse momentum interval (dpx) and the rapidity 
interval (dy) with all corrections applied. 

In Fig. [BJ the red dots and blue triangles represent the 
results from Arml and Arm2, respectively. The error 



bars and shaded rectangles indicate the one standard de- 
viation statistical and total systematic uncertainties, re- 
spectively. The total systematic uncertainties are given 
by adding all uncertainty terms except for the luminos- 
ity in quadrature. The vertical dashed lines shown in 
the rapidity range below 9.2 indicate the pt threshold of 
the Arm2 detector owing to the photon energy threshold 
and the geometrical acceptance. The px threshold of the 
Arml detector occurs at a higher value of pt than Arm2 
due to its smaller acceptance. A general agreement be- 
tween the Arml and Arm2 px spectra within statistical 
and systematic uncertainties is evident in Fig. [5] 

Fig. [7] presents the combined px spectra of the Arml 
and Arm2 detectors (black dots). The 68% confidence 
intervals incorporating the statistical and systematic un- 
certainties are indicated by the shaded green rectangles. 
The combined spectra below the pt threshold of Arml 
are taken from the Arm2 spectra alone. Above the px 
threshold of Arml, experimental px spectra of the Arml 
and Arm2 detectors have been combined following the 
"pull method" [44| and the combined spectra have ac- 
cordingly been obtained by minimizing the value of the 
chi-square function defined as 

X = 2^ — + ^penalty, (?) 

2—1 a— 1 \ / 

where the index i represents the px bin number running 
from 1 to n (the total number of px bins), N°\ s is the 
number of events and a a ,i is the uncertainty of the Arm-a 
analysis calculated by quadratically adding the statistical 
uncertainty and the energy scale uncertainty estimated 
by test beam data at SPS. The S a ,i denotes the system- 
atic correction to the number of events in the i-th bin of 
Arm-a: 

6 

s a ,i=Y,fiA- ( 8 ) 

The coefficient f 3 a i is the systematic shift of z-th bin con- 
tent due to the j-th systematic uncertainty term. The 
systematic uncertainty is assumed fully uncorrelated be- 
tween the Arml and Arm2 detectors, and consists of six 
uncertainties related to energy scale owing to the invari- 
ant mass shift, PID, beam center position, single-hit, po- 
sition dependent correction, and contamination by mulit- 
hit 7T° events. Coefficients e 3 a , which should follow a 
Gaussian distribution, can be varied to achieve the min- 
imum x 2 value in each chi-square test, while they are 
constrained by the penalty term 

6 

Xpcnalty = (l e Arml| 2 + ^Arn^l") ' (9) 

3=1 

The 7T° production rates for the combined data of LHCf 
are summarized in Tables IPV1 - IIX1 Note that the uncer- 
tainty in the luminosity determination ±6.1%, that is not 
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FIG. 6: (color online). Experimental pr spectra of the Arml (red dots) and Arm2 (blue triangles) detector. Error bars indicate 
the statistical uncertainties and shaded rectangles show the systematic uncertainties of the Arml and Arm2 detectors. 



included in Fig.[7J can make a pt independent shift of all 
spectra. 

For comparison, the px spectra predicted by various 
hadronic interaction models are also shown in Fig. [7] 
The hadronic interaction models that have been used 
in Fig. [7J arc dpmjet 3.04 (solid, red), QGSJET 11-03 
(dashed, blue), SIBYLL 2.1 (dotted, green), EPOS 1.99 
(dashed dotted, magenta), and pythia 8.145 (default 
parameter set, dashed double-dotted, brown). In these 
MC simulations, 7T°s from short lived particles that de- 
cay within 1 m from IP1, for example r\ — > 37T°, are also 
counted to be consistent with the treatment of the exper- 
imental data. Note that, since the experimental px spec- 
tra have been corrected for the influences of the detector 
responses, event selection efficiencies and geometrical ac- 
ceptance efficiencies, the px spectra of the interaction 
models may be compared directly to the experimental 
spectra as presented in Fig. [JJ 

Fig.[8]presents the ratios of px spectra predicted by the 
various hadronic interaction models to the combined px 
spectra. Error bars have been taken from the statistical 
and systematic uncertainties. A slight step found around 
Pt = 0.3 GeV in 8.9 < y < 9.0 is due to low pt cutoff of 
the Arml data. The ratios are summarized in Tables IXl- 



VII. DISCUSSION 

A. Transverse momentum spectra 

Several points can be made about Fig. [8] First, dpm- 
jet 3.04 and pythia 8.145 show overall agreement with 
the LHCf data for 9.2 < y < 9.6 and p T < 0.2 GeV, 
while the expected 7T° production rates by both models 
exceed the LHCf data as px becomes large. The lat- 
ter observation can be explained by the baryon/meson 
production mechanism that has been employed in both 
models. More specifically, the "popcorn model" [H, |46| 
is used to produce baryons and mesons through string 
breaking, and this mechanism tends to lead to hard pion 
spectra. SIBYLL 2.1, which is also based on the popcorn 
model, also predicts harder pion spectra than the exper- 
imental data, although the expected 7r° yield is generally 
small. 

On the other hand, QGSJET 11-03 predicts n° spec- 
tra that are softer than the LHCf data and the other 
models. This might be due to the fact that only one 
quark exchange is allowed in the QGSJET model. The 
remnants produced in a proton-proton collision are like- 
wise baryons with relatively small mass, so fewer pions 
with large energy are produced. 

Among hadronic interaction models tested in this anal- 




FIG. 7: (color online). Combined pt spectra of the Arml and Arm2 detectors (black dots) and the total uncertainties (shaded 
rectangles) compared with the predicted spectra by hadronic interaction models. 



ysis, EPOS 1.99 shows the best overall agreement with the 
LHCf data. However EPOS 1.99 behaves softer than the 
data in the low px region, px ^ 0.4 GeV in 9.0 < y < 9.4 
and pt < 0.3 GeV in 9.4 < y < 9.6, and behaves harder 
in the large px region. Specifically a dip found in the 
ratio of EPOS 1.99 to the LHCf data for y > 9.0 can be 
attributed to the transition between two pion produc- 
tion mechanisms: string fragmentation via cut Pomeron 
process (low energy ~ low px for the fixed rapidity) and 
remnants of projectile/target (high energy ~ large px for 
the fixed rapidity) [471 ]. 



B. Average transverse momentum 

According to the scaling law proposed by several au- 
thors [4a - l50j . the average transverse momentum as a 
function of rapidity should be independent of the center 
of mass energy in the projectile fragmentation region. 
Average transverse momentum, (px), can be obtained 
by fitting an empirical function to the px spectra in each 
rapidity range. In this analysis, among several ansatz 
proposed for fitting the px spectra, an exponential dis- 
tribution has been first chosen with the form 



Omcl dp 3 



This distribution is motivated by a thermodynamical 
model [HTJ] . The parameter A [GeV -2 ] is a normaliza- 
tion factor and T [GeV] is the temperature of 7r°s with 
a given transverse momentum px. Using Eq. 1101 (px) is 
derived as a function of T: 



(pi 



2 K 3/2 (m^/T) : 



111) 



^■exp(-Jp T 2 +m 2 /T). 



(10) 



where K a (m T o /T) is the modified Bessel function. 

Best-fit results for T and (px) are summarized in 
Table [TTJ The worse fit quality values are found for 
9.2 < y < 9.4 (r7 d of = 3.6) and 9.4 < y < 9.6 
(% 2 /dof = 11.1). These are caused by data points near 
Pt = 0.25 GeV which exceed the best-fit exponential 
distribution and the experimental px spectra decreas- 
ing more rapidly than Eq. (JTOJl for p T > 0.3 GeV. The 
upper panels in Fig. |H] show the experimental px spec- 
tra (black dots and green shaded rectangles) and the 
best-fit of Eq. (fl"U|) (dashed curve) in the rapidity range 
9.2 < y < 9.4 and 9.4 < y < 9.6. The bottom panels in 
Fig. IH1 show the ratio of the best-fit distribution to the 
experimental data (blue triangles). Shaded rectangles in- 
dicate the statistical and systematic uncertainties. Even 
though the minimum x 2 /dof values are large, the best- 
fit T values are consistent with temperatures that are 
typical of soft QCD processes and the predictions of the 
thermodynamical model (T < 180 MeV) [U for y > 8.9. 
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FIG. 8: (color online). Ratio of the combined pt spectra of the Arml and Arm2 detectors to the predicted pt spectra by 
hadronic interaction models. Shaded areas indicate the range of total uncertainties of the combined pt spectra. 



Another possibility is that thepx distributions in Fig. [7] 
can also be described by a Gaussian distribution: 



CTincl dp 3 



(12) 



Gauss 



The Gaussian width ercauss determines the mean square 
Pt of the pt spectra, (px) is derived as a function of 
CTGauss according to: 



(pt) 



j2p T f(p T )dp T = VtF_ 
/ 2p T f(p T )dp T 2 



-^Ga 



(13) 



where /(px) is given by Eq. (|12l) . Best-fit results for 
CGauss and (px) are summarized in Table [TXJ In this case 
good fit quality values are found for all rapidity ranges. 
The best-fit of Eq. (fT2"j) (dotted curve) and the ratio of the 
best-fit Gaussian distribution to the experimental data 
(red open boxes) are found in Fig. [9] 

A third approach for estimating (px) is simply numeri- 
cally integrating the px spectra. With this approach (px) 
is given by 



(Pt) 



J °° 2Trp T f(p T )dp T 
/ oo 27rp T /(px)dpx 



(14) 



where /(px) is the measured spectrum given in Fig. [7] for 
each of the six ranges of rapidity. In this analysis, (px) 



is obtained over the rapidity range 9.2 < y < 11.0 where 
the pt spectra are available down to GeV. Although the 
upper limits of numerical integration are actually finite, 
^upper ^ 0.6 GeV, the contribution of the high pt tail 
to (pt) is negligible. p T ppcl and the obtained (px) are 
summarized in Table [TT1 



The values of (px) obtained by the three methods dis- 
cussed above are in general agreement. When a specific 
values of (pt) are needed for this paper the values chosen 
((.PT)LHCf) are defined as follows. For the rapidity range 
8.9 < y < 9.2, (px)LHCf is taken from the weighted mean 
of (pt) obtained by the exponential fit of Eq. (|lll) and 
the Gaussian fit of Eq. (fl"3|) . The systematic uncertainty 
related to a possible bias of the (px) extraction methods 
is estimated by the difference of (pt) derived from these 
two different fitting functions. The estimated systematic 
uncertainty is ±6 % for both rapidity bins. For the ra- 
pidity range 9.2 < y < 11.0, the results obtained by the 
Gaussian fit and numerical integration are used to cal- 
culate the weighted mean of (px)LHCf in order to avoid 
the poor quality of fit of the exponential function in this 
rapidity range. Systematic uncertainty is estimated to be 
±3 % and ±2 % for 9.2 < y < 9.4 and 9.4 < y < 11.0, re- 
spectively. The values of (px)LHCf obtained by the above 
calculation are summarized in Table IIIII 
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FIG. 9: (color online). (Upper) Experimental px spectra 
(black dots and green shaded rectangles), the best-fit expo- 
nential distributions (Eq. (|10[) , dashed curve) and the best-fit 
Gaussian distributions (Eq. 1)120 . dotted curve). (Bottom) 
Ratios of the best-fit exponential or Gaussian distribution to 
the experimental data (blue triangles or red open boxes) and 
the statistical and systematic uncertainties (green shaded ar- 
eas). For both the upper and bottom panels, the rapidity 
ranges 9.2 < y < 9.4 and 9.4 < y < 9.6 are shown on the left 
and right panels, respectively. 





Exponential fit 


Gaussian fit 


Numerical integration 


Rapidity 


X* (dof) 


T (pr> Stat, error 


X* (dof) 


""Gauss (Pt) 


Stat, error 


upper 

Pt 


<Pt) 


Stat, error 






[MeV] [MeV] 


[MeV] 




[MeV] [MeV] 


[MeV] 


[GeV] [MeV] 


[MeV] 


[8.9, 9.0] 


0.6 (7) 


83.8 201.4 


13.5 


2.0 (7) 


259.0 229.6 


13.1 








[9.0, 9.2] 


8.2 (7) 


75.2 184.1 


5.0 


0.9 (7) 


234.7 208.0 


4.6 








[9.2, 9.4] 


28.7 (8) 


61.7 164.0 


2.8 


6.9 (8) 


201.8 178.9 


3.4 


0.6 


167.7 


9.6 


[9.4, 9.6] 


66.3 (6) 


52.8 140.3 


1.9 


3.3 (6) 


166.3 147.4 


2.7 


0.4 


144.8 


3.2 


[9.6, 10.0] 


14.0 (5) 


43.3 123.5 


2.2 


0.3 (5) 


139.2 123.3 


3.0 


0.4 


117.0 


2.1 


[10.0, 11.0] 


9.0 (2) 


21.3 77.7 


2.3 


2.1 (2) 


84.8 75.1 


2.9 


0.2 


76.9 


2.6 



TABLE II: Best-fit results of exponential and Gaussian px functions to the LHCf data and average n transverse momenta for 
the rapidity range 8.9<y<11.0 obtained by using the exponential fit, Gaussian fit and numerical integration. 



The values of (pt) that have been obtained in this 
analysis, shown in Table IIII1 are compared in Fig. [TO] 
with the results from UA7 at SppS (y/s = 630 GeV) @ 
and the predictions of several hadronic interaction mod- 
els. In Fig. [TO] (px) is presented as a function of rapidity 
loss Ay = 2/beam - U, where beam rapidity y b cam is 8.92 
for Vs = 7TcV and 6.50 for y/s = 630 GeV. This shift 
of rapidity scales the results with beam energy and it al- 
lows a direct comparison between LHCf results and past 
experimental results at different collision energies. The 



black dots and the red diamonds indicate the LHCf data 
and the UA7 results, respectively. Although the LHCf 
and UA7 data in Fig. [TO] have limited overlap and the 
systematic errors of the UA7 data are relatively large, 
the (pr) spectra for LHCf and UA7 in Fig. [TO] mostly 
appear to lie along a common curve. 

The (pr) predicted by hadronic interaction models are 
shown by open circles (sibyll 2.1), open boxes (QGSJET 
11-03) and open triangles (epos 1.99). SIBYLL 2.1 typi- 
cally gives harder w° spectra (larger (pt)) and QGSJET 
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Rapidity 


(pt) 


Total uncertainty 




[MeV] 


[MeV] 


[8.9, 9.0] 


215.3 


17.3 


[9.0, 9.2] 


196.8 


12.5 


[9.2, 9.4] 


172.2 


5.9 


[9.4, 9.6] 


146.3 


3.9 


[9.6, 10.0] 


119.2 


3.4 


[10.0, 11.0] 


75.8 


2.9 



TABLE III: Average transverse momentum of n for the ra- 
pidity range 8.9<y<11.0. Total pt uncertainty includes both 
the statistical and systematic uncertainties. 



11-03 gives softer ir° spectra (smaller (pt)) than the ex- 
perimental data. For each prediction, solid and dashed 
lines indicate (pt) &t the center of mass energy at SppS 
and the LHC, respectively. Of the three models the pre- 
dictions by EPOS 1.99 show the smallest dependence of 
(px) on the two center of mass energies, and this tendency 
is consistent with the LHCf and UA7 results except for 
the UA7 data at Ay = —0.15 and 0.25. It is also evi- 
dent in Fig. [TU] that amongst the three models the best 
agreement with the LHCf data is obtained by EPOS 1.99. 




FIG. 10: (color online). Average pt as a function of rapid- 
ity loss Ay. Black dots and red diamonds indicate the LHCf 
data and UA7 results taken from Ref. [f|, respectively. The 
predictions of hadronic interaction models are shown by open 
boxes (sibyll 2.1), open circles (qgsjet 11-03) and open tri- 
angles (epos 1.99). For the predictions of the three models, 
solid and dashed curves indicate the results for the center of 
mass energy at the SppS and the LHC, respectively. 



VIII. CONCLUSIONS 

The inclusive production of neutral pions in the rapid- 
ity range larger than y = 8.9 has been measured by the 
LHCf experiment in proton-proton collisions at the LHC 
in early 2010. Transverse momentum spectra of neu- 
tral pions have been measured by two independent LHCf 
detectors, Arml and Arm2, and give consistent results. 
The combined Arml and Arm2 spectra have been com- 
pared with the predictions of several hadronic interaction 
models, dpmjet 3.04, epos 1.99 and pythia 8.145 agree 
with the LHCf combined results in general for the rapid- 
ity range 9.0 < y < 9.6 and px < 0.2 GeV. QGSJET 11-03 
has poor agreement with LHCf data for 8.9 < y < 9.4, 
while it agrees with LHCf data for y > 9.4. Among the 
hadronic interaction models tested in this paper, epos 
1.99 shows the best overall agreement with the LHCf data 
even for y > 9.6. 

The average transverse momentum, (pt), of the com- 
bined px spectra is consistent with typical values for soft 
QCD processes. The (pt) spectra for LHCf and UA7 in 
Fig. [TU] mostly appear to lie along a common curve. The 
(Pt) spectra derived by LHCf agrees with the expectation 
of EPOS 1.99. Additional experimental data are needed 
to establish the dependence, or independence, of (pt) on 
the center of mass collision energy. 
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Appendix 

The inclusive production rates of 7r°s measured by 
LHCf arc summarized in Tables IIVI IIXI The ratios of 
inclusive production rates of 7r°s predicted by MC sim- 
ulations to the LHCf measurements are summarized in 
Tables HIXVl 
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Pt range 


Production rate 


Stat, uncertainty 


Syst.+Stat. uncertainty 


[GeV] 


[GeV~ 2 ] 


[GeV" 2 ] 


[GeV" 2 ] 


[0.10, 0.15] 


2. 71x10^ 


±1.41xl0 _i 


-1.41xl0-\ +1.58xl0 _i 


[0.15, 0.20] 


1.95X10" 1 


±8.85xl0- 2 


-8.85xl0- 2 , +9.95xl0- 2 


[0.20, 0.25] 


1.25xl0 _1 


±4.98xl0~ 2 


-4.98xl0" 2 , +5.66xl0" 2 


[0.25, 0.30] 


7.15xl0" 2 


±2.54xl0- 2 


-2.54xl0- 2 , +2.90xl0- 2 


[0.30, 0.35] 


4.34xl0" 2 


±3.21xl0~ 3 


-4.21X10- 3 , +4.22xl0" 3 


[0.35, 0.40] 


2.36xl0" 2 


±2.45xl0- 3 


-3.65xl0- 3 , +3.66xl0" 3 


[0.40, 0.45] 


1.50xl0" 2 


±2.05xl0- 3 


-3.25xl0" 3 , +3.26xl0" 3 


[0.45, 0.50] 


6.73xl0^ 3 


±1.50xl0~ 3 


-2.15xl0- 3 , +2.16xl0" 3 


[0.50, 0.60] 


3.42xl0" 3 


±8.08xl0" 4 


-1.27xl0" 3 , +1.27xl0" 3 


TABLE IV: 


Production rate for the 7r° 


production in the rapidity 


range 8.9<y<9.0. 



Pt range 


Production rate 


Stat, uncertainty 


Syst.+Stat. uncertainty 


[GeV] 


[GeV- 2 ] 


[GeV" 2 ] 


[GeV" 2 ] 


[0.10, 0.15] 


2.30X10- 1 


±8.06xl0" 2 


-8.06xl0- 2 , +8.12xl0~ 2 


[0.15, 0.20] 


1.51X10- 1 


±3.99xl0- 2 


-3.99xl0- 2 , +4.19xl0- 2 


[0.20, 0.25] 


8.92xl0" 2 


±2.62xl0- 3 


-3.14xl0- 3 , +3.14xl0" 3 


[0.25, 0.30] 


5.43xl0 -2 


±2.18xl0- 3 


-3.35xl0" 3 , +3.35xl0" 3 


[0.30, 0.35] 


3.21xl0- 2 


±1.79xl0- 3 


-3.31 xlO -3 , +3.32xl0- 3 


[0.35, 0.40] 


1.80xl0 -2 


±1.44xl0- 3 


-2.72xl0" 3 , +2.73xl0 -3 


[0.40, 0.45] 


8.91xl0- 3 


±1.00xl0- 3 


-1.83xl0- 3 , +1.84xl0- 3 


[0.45, 0.50] 


3.59xl0" 3 


±6.49xl0- 4 


-l.OlxlO- 3 , +1.01xl0 -3 


[0.50, 0.60] 


9.72xl0 -4 


±2.65 xlO -4 


-3.66xl0- 4 , +3.65xl0" 4 


TABLE V: 


Production rate for the 


7T° production in the rapidity 


range 9.0<y<9.2. 



Pt range 


Production rate 


Stat, uncertainty 


Syst.+Stat. uncertainty 


[GeV] 


[GeV" 2 ] 


[GeV" 2 ] 


[GeV" 2 ] 


[0.00, 0.05] 


3. 31x10-' 


±1.58xl0 _i 


-1.58xl0-\ +1.84X10- 1 


[0.05, 0.10] 


2.31xl0 _1 


±9.05xl0- 2 


-9.05xl0" 2 , +1.07X10" 1 


[0.10, 0.15] 


1.66X10- 1 


±3.23xl0- 3 


-1.96xl0- 2 , +1.96xl0- 2 


[0.15, 0.20] 


1.06xl0 _1 


±2.42xl0- 3 


-3.76xl0" 3 , +3.76xlO" 3 


[0.20, 0.25] 


5.71xl0- 2 


±1.91xl0- 3 


-2.69xl0- 3 , +2.69xl0- 3 


[0.25, 0.30] 


3.58xl0" 2 


±1.65xl0- 3 


-2.97xl0" 3 , +2.98xl0 -3 


[0.30, 0.35] 


1.77xl0- 2 


±1.26xl0- 3 


-2.29xl0- 3 , +2.29xl0" 3 


[0.35, 0.40] 


8.07xl0 -3 


±9.02xl0- 4 


-1.49xl0 -3 , +1.49xl0 -3 


[0.40, 0.50] 


1.35xl0- 3 


±2.66xl0- 4 


-3.71xl0- 4 , +3.72xl0- 4 


[0.50, 0.60] 


1.47xl0- 4 


±8.16xl0 -5 


-9.17xl0" 5 , +9.18xl0" 5 


TABLE VI: 


Production rate for the n 


production in the rapidity 


range 9.2<y<9.4. 



Pt range 


Production rate 


Stat, uncertainty 


Syst.+Stat. uncertainty 


[GeV] 


[GeV" 2 ] 


[GeV" 2 ] 


[GeV" 2 ] 


[0.00, 0.05] 


2. 03x10-' 


±8.63xl0" 3 


-3.09xl0- 2 , +3.09xl0" 2 


[0.05, 0.10] 


1.73X10- 1 


±3.75xl0- 3 


-1.64xl0- 2 , +1.64xl0- 2 


[0.10, 0.15] 


1.07xl0 _1 


±2.24xl0- 3 


-4.61xl0~ 3 , +4.60xl0 -3 


[0.15, 0.20] 


6.30xl0- 2 


±1.80xl0- 3 


-2.45xl0- 3 , +2.45xl0- 3 


[0.20, 0.25] 


3.20xl0 -2 


±1.51xl0- 3 


-2.63xl0 -3 , +2.64xl0 -3 


[0.25, 0.30] 


1.45xl0- 2 


±1.17xl0- 3 


-2.14xl0- 3 , +2.15xl0- 3 


[0.30, 0.35] 


3.64xl0 -3 


±6.44xl0 -4 


-9.28xl0 -4 , +9.29xl0 -4 


[0.35, 0.40] 


1.54xl0" 3 


±4.88xl0- 4 


-6.20xl0- 4 , +6.21xl0- 4 


[0.40, 0.50] 


5.43xl0 -5 


±6.19xl0 -5 


-6.41xl0" 5 , +6.41xl0" 5 



TABLE VII: Production rate for the ir° production in the rapidity range 9.4<y<9.6. 
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p T range 
[GeV] 


Production rate 
[GeV~ 2 ] 


Stat, uncertainty 
[GeV -2 ] 


Syst.+Stat. uncertainty 
[GeV- 2 ] 


[0.00, 0.05] 


1. 20x10^ 


±3.49xl0 -3 


-9.66xl0" 3 , +9.68xl0~ 3 


[0.05, 0.10] 


8.28xl0 -2 


±1.55xl0 -3 


-2.89X10 -3 , +2.90X10 -3 


[0.10, 0.15] 


4.49xl0" 2 


±1.02xl0~ 3 


-1.88xl0~ 3 , +1.88xl0" 3 


[0.15, 0.20] 


2.10xl0 -2 


±8.40xl0 -4 


-1.28X10 -3 , +1.28X10 -3 


[0.20, 0.25] 


7.43xl0" 3 


±6.05xl0~ 4 


-9.73xl0" 4 , +9.76xl0 -4 


[0.25, 0.30] 


1.84xl0 -3 


±4.05 xlO -4 


-5.15X10 -4 , +5.16xl0 -4 


[0.30, 0.40] 


2.17xl0" 4 


±1.21xl0" 4 


-1.33xl0" 4 , +1.33xl0" 4 



TABLE VIII: Production rate for the ir production in the rapidity range 9.6<y<10.0. 



Pt range 


Production rate 


Stat, uncertainty 


Syst.+Stat. uncertainty 


[GeV] 


[GeV" 2 ] 


[GeV" 2 ] 


[GeV" 2 ] 


[0.00, 0.05] 


1.28xl0- 2 


±9.69xl0 -4 


-1.42xl0- 3 , +1.42xl0" 3 


[0.05, 0.10] 


7.55xl0~ 3 


±3.79xl0- 4 


-8.88xl0- 4 , +8.85xl0- 4 


[0.10, 0.15] 


2.37xl0" 3 


±1.95xl0- 4 


-3.77xl0 -4 , +3.76xl0 -4 


[0.15, 0.20] 


1.91xl0 -4 


±6.22xl0- 5 


-6.99X10- 5 , +6.98xl0- 5 


[0.20, 0.30] 


8.37xl0 -6 


±1.03xl0" 5 


-1.05xl0" 5 , +1.05xlQ- 5 



TABLE IX: Production rate for the 7r° production in the rapidity range 10.0<y<11.0. 



Pt range 
[GeV] 


DPMJET 3.04 


QGSJET 11-03 


SIBYLL 2.1 


epos 1.99 


PYTHIA 8.145 


[0.10, 0.15] 


1.36 


1.37 


0.74 


1.23 


1.38 


[0.15, 0.20] 


1.59 


1.48 


0.85 


1.41 


1.57 


[0.20, 0.25] 


1.97 


1.71 


1.04 


1.79 


2.03 


[0.25, 0.30] 


1.82 


1.34 


1.00 


1.57 


2.02 


[0.30, 0.35] 


1.32 


0.71 


0.77 


1.04 


1.53 


[0.35, 0.40] 


1.57 


0.69 


0.97 


1.02 


1.87 


[0.40, 0.45] 


1.70 


0.56 


1.08 


0.83 


1.89 


[0.45, 0.50] 


2.54 


0.59 


1.60 


1.01 


2.81 


[0.50, 0.60] 


2.76 


0.38 


1.73 


0.90 


3.05 



TABLE X: Ratio of 7r production rate of MC simulation to data in the rapidity range 8.9<y<9.0. 



PT range 
[GeV] 


DPMJET 3.04 


QGSJET 11-03 


SIBYLL 2.1 


epos 1.99 


PYTHIA 8.145 


[0.10, 0.15] 


1.23 


1.20 


0.64 


1.13 


1.26 


[0.15, 0.20] 


1.23 


1.06 


0.62 


1.09 


1.30 


[0.20, 0.25] 


1.11 


0.81 


0.60 


0.88 


1.28 


[0.25, 0.30] 


1.14 


0.68 


0.64 


0.94 


1.34 


[0.30, 0.35] 


1.27 


0.58 


0.73 


0.88 


1.44 


[0.35, 0.40] 


1.51 


0.52 


0.84 


0.76 


1.65 


[0.40, 0.45] 


2.08 


0.47 


1.08 


0.74 


2.15 


[0.45, 0.50] 


3.43 


0.53 


1.69 


0.93 


3.33 


[0.50, 0.60] 


6.43 


0.48 


2.75 


1.45 


5.82 



TABLE XI: Ratio of ir production rate of MC simulation to data in the rapidity range 9.0<y<9.2. 
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Pt range 
[GeV] 


DPMJET 3.04 


QGSJET 11-03 


SIBYLL 2.1 


epos 1.99 


PYTHIA 8.145 


[0.00, 0.05] 


1.07 


1.29 


0.54 


0.98 


1.04 


[0.05, 0.10] 


1.22 


1.33 


0.60 


1.10 


1.23 


[0.10, 0.15] 


1.02 


0.95 


0.48 


0.89 


1.07 


[0.15, 0.20] 


0.97 


0.78 


0.48 


0.77 


1.07 


[0.20, 0.25] 


1.10 


0.70 


0.55 


0.86 


1.24 


[0.25, 0.30] 


1.12 


0.52 


0.56 


0.78 


1.21 


[0.30, 0.35] 


1.50 


0.47 


0.70 


0.66 


1.48 


[0.35, 0.40] 


2.20 


0.41 


0.83 


0.66 


1.93 


[0.40, 0.50] 


6.76 


0.59 


1.99 


1.46 


5.37 


[0.50, 0.60] 


15.90 


0.26 


3.34 


3.12 


11.36 



TABLE XII: Ratio of n production rate of MC simulation to data in the rapidity range 9.2<y<9.4. 



Pt range 
[GeV] 


DPMJET 3.04 


QGSJET 11-03 


SIBYLL 2.1 


epos 1.99 


PYTHIA 8.145 


[0.00, 0.05] 


1.11 


1.28 


0.51 


0.96 


1.14 


[0.05, 0.10] 


0.97 


1.00 


0.44 


0.84 


1.00 


[0.10, 0.15] 


1.00 


0.89 


0.46 


0.77 


1.07 


[0.15, 0.20] 


1.02 


0.71 


0.46 


0.76 


1.14 


[0.20, 0.25] 


1.27 


0.63 


0.55 


0.84 


1.28 


[0.25, 0.30] 


1.82 


0.54 


0.66 


0.72 


1.60 


[0.30, 0.35] 


4.71 


0.74 


1.32 


1.09 


3.68 


[0.35, 0.40] 


6.60 


0.48 


1.39 


1.28 


4.79 



TABLE XIII: Ratio of 7r production rate of MC simulation to data in the rapidity range 9.4<y<9.6. 



Pt range 
[GeV] 


DPMJET 3.04 


QGSJET 11-03 


SIBYLL 2.1 


epos 1.99 


PYTHIA 8.145 


[0.00, 0.05] 


0.98 


1.07 


0.39 


0.75 


1.00 


[0.05, 0.10] 


1.04 


0.96 


0.42 


0.79 


1.09 


[0.10, 0.15] 


1.23 


0.86 


0.49 


0.88 


1.24 


[0.15, 0.20] 


1.66 


0.73 


0.56 


0.87 


1.45 


[0.20, 0.25] 


2.96 


0.67 


0.71 


0.87 


2.23 


[0.25, 0.30] 


6.42 


0.72 


1.07 


1.26 


4.30 


[0.30, 0.40] 


14.86 


0.65 


1.51 


2.27 


8.85 



TABLE XIV: Ratio of tt production rate of MC simulation to data in the rapidity range 9.6<y<10.0. 



Pt range 


DPMJET 3.04 


QGSJET 11-03 


SIBYLL 2.1 


epos 1.99 


PYTHIA 8.145 


[GeV] 












[0.00, 0.05] 


2.12 


1.15 


0.52 


1.03 


1.66 


[0.05, 0.10] 


2.46 


1.00 


0.54 


0.98 


1.79 


[0.10, 0.15] 


4.25 


1.05 


0.67 


1.09 


2.80 


[0.15, 0.20] 


22.24 


2.47 


2.22 


3.41 


13.42 



TABLE XV: Ratio of tt° production rate of MC simulation to data in the rapidity range 10.0<y<11.0. 
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